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Abstract A strong coupling between convection and 
pulsations is known to play a major role in the dis- 
appearance of unstable modes close to the red edge of 
the classical Cepheid instability strip. As mean-field 
models of time-dependent convection rely on weakly- 
constrained parameters (e.g. Baker 1987), we tackle this 
problem by the means of 2-D Direct Numerical Simu- 
lations (DNS) of A>:-mechanism with convection. 

Using a linear stability analysis, we first determine 
the physical conditions favourable to the /^-mechanism 
to occur inside a purely-radiative layer. Both the in- 
stability strips and the nonlinear saturation of unstable 
modes are then confirmed by the corresponding DNS. 
We next present the new simulations with convection, 
where a convective zone and the driving region over- 
lap. The coupling between the convective motions and 
acoustic modes is then addressed by using projections 
onto an acoustic subspace. 

Keywords Hydrodynamics - Instabilities - Stars: os- 
cillations - Convection - Methods: numerical 



1 Introduction 



EddingtonI (|l917l ) discovered an excitation mechanism 
of stellar oscillations that is related to the opacity be- 
haviour in ionisation zones: the mechanism, where 
K denotes the opacity. This mechanism occurs when 
the opacity varies during compressio n phases so as 
to block the emerg:ing: radiative flux ( ZhevakinI 1963 : 
Cox and Whitnev 19581 ) . Ionisation regions correspond 
to a strong increase in opacity that leads to the opac- 
ity bumps responsible for the local excitation of modes. 
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These driving zones should nevertheless be located at 
a precise radius inside the star, neither too close to 
the surface nor to deep, in order to balance the damp- 
ing that mainly occurs at the surface. It defines the 
so-called transition region that shapes the limit be- 
tween the quasi-adiabatic interior and the strongly non- 
adiabatic surface. For classical Cepheids that pul- 
sate in the fundamental acoustic mode, this transi- 
tion region is located at a temperature T 4 x 
10^ K corresponding to the second helium ionisation 
([Baker and Kippenhahnlll965l ). However, the bump lo- 
cation is not solely responsible for the acoustic instabil- 
ity. A careful treatment of the A>:-mechanism would in- 
volve dynamical couplings with convection, metallicity 
effects and both realistic equations of state and opac- 
ity tables ( Bono, Marconi, and Stellingwert 1999). The 
purpose of our model is to simplify the hydrodynamic 
approach while retaining the leading- order phenomenon 
-e.g. the location of the opacity bump or its amplitude- 
such that feasible direct numerical simulations of the 
mechanism with convection can be achieved. 



2 The first step: radiative models in 1-D 

We first focus on radial modes propagating in a purely- 
radiative and partially-ionised layer and then restrict 
our study to the 1-D case. Our model consists of a 
local zoom about an ionisation region and is composed 
by a monatomic and perfect gas (with 7 = Cp/cy =5/3) 
under both a constant gravity g and kinematic viscosity 
u. All quantities in our model are dimensionless, e.g. 
temperature is given in unit of the surface temperature 
of the star, length in unit of a fraction of the star radius 
and velocity in unit of the sound speed at the surface. 

In order to easily investigate the influence of the 
opacity bump on the stability, we adopt the following 
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Fig. 1 Influence of the hollow parameters on the conduc- 
tivity profile for Kmax = 10"^ and Tbump/^surf = 3.5. 



co nductivity hollow to mimic th is bump (as hi oc 1/i^, 
see 



Gastine and Dintrans l2QQ8al . hereafter Paper I): 



-7r/2 + arctan(crT+T-) 



7r/2 + arctan(cre^) 



, (1) 



with A (i^max-^min)/^max and - Tbump ± 

e. Here Tq is the equilibrium temperature profile, Tbump 
is the hollow central temperature, and <j, e and A de- 
note its slope, width and relative amplitude, respec- 
tively. Examples of common values of these parameters 
are provided in Fig. [T] 

2.1 Conditions for the instability 

Using the well-k nown work integral formalism (e.g. 
Unno et a/.l 119891 ). it is easy to show that the follow- 
ing condition is necessary to get unstable modes by hz- 
mechanism in variable stars: 



dlCr 
dz 



< where JCi 



dlnKp 
d\nTo 



(2) 



However, this condition is not sufficient as the ionisa- 
tion region should also be located at a given location 
in the star, neither too deep nor too close to its sur- 
face. This area is the so-called transition region that 
separates the quasi-adiabatic interior from t he str ongly 
non-adiabatic surface and is defined by (e.g. Ic^lTisQi ) 



{cyTo)Am 
UL 



(3) 



where Am is the integrated mass between the consid- 
ered point and the surface, H the mode period and L the 
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Fig. 2 a) Vertical profile of the ^ coefficient (Eq. [S). 
The superimposed green zone represents the location where 
^ = 1 ± 0.5. b) The real part of the work integral plotted 
for the three different equilibrium models. The dot-dashed 
black line is for the constant radiative conductivity case, 
c) Corresponding radiative conductivity profiles, d) Corre- 
sponding equilibrium fields dJCr/dz (Eq. |2l). 



luminosity. ^ is the ratio between the thermal energy 
embedded between the given radius and the surface and 
the energy radiated during an oscillation period. 

2.2 The linear stability analysis 

In order to check these two conditions, we perform a lin- 
ear stabil ity analysis b y using the Linear Solver Builder 
code ( V aldettaro et al. 2007). By starting with the con- 
ductivity profiles shown in Fig. [TJ we first compute the 
equilibrium setup that is solution of both the hydro- 
static and radiative equilibria given by 



Vpo = P^g and div[i^o(^o)VTo] = 0, 



(4) 



where po and po denote the equilibrium pressure and 
density. Then, we solve the linear oscillation equations 
for the perturbations by seeking normal modes of the 
form exp(At) with A = r + zcj, where uo is the mode 
frequency and r its growth or damping rate (unstable 
modes corresponding to r > 0). By varying the hol- 
low parameters, we then determine the physical con- 
ditions that lead to unstable modes excited by the k- 
mechanism in our layer (cf. Paper I). 

The two criteria (0 and ([3]) predict that modes are 
unstable for a "sufficient" hollow (adequate amplitude 
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and width) that is located in the transition region. 
Using the work integral computed from the obtained 
eigenvectors, we indeed recover instability strips for the 
Ai:-mechanism by considering the three different equilib- 
rium setups (Fig. [2|): 

• The first case (dotted blue line, Tbump/^surf = 1-7) 
corresponds to a hot star with an ionisation region 
close to the surface. As the surrounding density is 
small, ^ <C 1 and the conductivity hollow hardly in- 
fluences the work integral: driving is unable to prevail 
over damping, leading to r < or stable modes. 

• In the second case (solid green line, Tbump/^surf = 
2.1), the radiative conductivity begins to decrease 
significantly at the location of the transition region 
where ^ 1. As a consequence, driving becomes im- 
portant. Furthermore, the radiative conductivity in- 
crease occurs where non-adiabatic effects are already 
significant, i.e. ^ < 1 there. It means that no damp- 
ing occurs between the hollow position and the sur- 
face as the radiative flux perturbations are frozen-in. 
Driving is overcoming damping in this case, therefore 
r > and modes are unstable. 

• The third case (dashed red line, Tbump/^surf = 2.8) 
corresponds to a cold star of which ionisation re- 
gion is located deeper in the stellar atmosphere where 
^ ^ 1. Ionisation then occurs in a quasi- adiabatic 
location. As a consequence, the excitation provided 
by the conductivity hollow cannot balance the damp- 
ing arising at the top of the layer, thus r < and 
modes are stable. 

In conclusion, our first radiative model of the k,- 
mechanism is well able to reproduce the main physics 
involved in the oscillations of Cepheid stars: (i) one 
should have dJCr/dz < to drive the oscillations; (ii) 
the thermal engine underlying the conductivity hollow 
should also be located in the transition region where 
^ 1 (Fig.[2|). 

2.3 Study of the nonlinear saturation 

To confirm both the results obtained previously in the 
linear stability analysis and to investigate the result- 
ing nonlinear saturation, we perform direct numerical 
simulations of our problem. The idea is to start from 
the most favourable initial conditions found in the sta- 
bility analysis and then advance the general hydro- 
dynamic equations in time. We use a modified ver- 
sion of the public-domain finite-difference Pencil Cod^ 
that includes an implicit solver of our own for the 
radiative diffusion term in t he temperature equation 
( Gastine and Dintrans 2QQ8bl . hereafter Paper II). 




^See ,http : //pencil- code . googlecode . com/| 
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Fig. 3 a) Temporal power spectrum for the momentum 
in the {z, C(;)-plane. b) The resulting spectrum after in- 
tegrating over depth, c) Comparison between normalised 
momentum profiles for n = (0, 2) modes according to the 
DNS power spectrum (solid black line) and to the linear 
stability analysis (dotted blue line). 



To determine which modes are present in the DNS 
in the nonlinear-saturation regime, we first perform 
a temporal Fourier transform of the momentum field 
pu{z^ t) and plot the resulting power spectrum in the 
(z,a;)-plane (Fig. |3ti). The mean spectrum follows by 
integrating pu{z^uj) over depth (Fig. ISJd). With this 
method, acoustic modes are extracted because they 
emerge a s ^'shark-fin profiles" about definite eigenfre- 
quencies ( Dintrans and Brandenbur3 [2QQ4). 

Several discrete peaks corresponding to normal 
modes well appear but the fundamental mode close 
to cuq = 5.439 clearly dominates. Moreover, the lin- 
ear eigenfunctions are compared to the mean profiles 
computed from a zoom taken in the DNS power spec- 
trum about eigenfrequencies uq = 5.439 and UJ2 = 11.06 
(Fig. [3t). The agreement between the linear eigenfunc- 
tions (dotted blue lines) and the DNS profiles (solid 
black lines) is remarkable. In summary. Fig. |3] shows 
that several overtones are present in the DNS, even for 
long times. These overtones are however linearly stable 
suggesting that some underlying energy transfers occur 
between modes through nonlinear couplings. 

These nonlinear interactions are investigated thanks 
to a tool that nieasures the sound generation by turbu - 
lent flows (e.g. Bogdan. Cattaneo. and Malagolil 19931 ). 
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It involves projections of the DNS physical fields onto 
the regular and adjoint eigenvectors that are solutions 
to the linear oscillation equations. The projection co- 
efficients then give the time evolution of each acoustic 
mode that is present in the DNS separately. Further- 
more, as the evolution of the kinetic energy of each 
mode is also accessible with this method, one can fol- 
low the energy transfer between modes. 




Fig. 4 Time evolution of the kinetic energy ratio for the 
n = and n = 2 acoustic modes. 

Figure [4] shows an example of the time evolution 
of the kinetic energy ratio En/E^^ for the two modes 
n = (i.e. the fundamental one) and n = 2 (i.e. the 
second overtone). Here E^ denotes the kinetic energy 
embedded in the acoustic mode of order n, while E^^ is 
the total kinetic energy in the simulation. After its lin- 
ear transient growth, a given fraction of the fundamen- 
tal mode energy is progressively transferred to the sec- 
ond overtone until the nonlinear saturation is achieved 
above time t o:^ 150, with finally 85% of the total kinetic 
energy in the n = mode and the remaining 15% in the 
n = 2 one. We recall that it is at a first glance surpris- 
ing to detect this second overtone at long times because 
it is linearly stable and should normally be damped by 
diffusive effects. 

The reason for its presence lies in a favored coupling 
in the period ratio between these two modes. Indeed, 
the fundamental period is Pq = 27r/cJo — 1.155 in this 
simulation while the n = 2 one is P2 — 0.568. As a 
consequence, the corresponding period ratio is close to 
one half {P2/P0 — 0.491) such that the second over- 
tone is involved in the nonlinear saturation through a 
2:1 resonance with the fundamental mode. The lin- 
ear growth of the fundamental mode is then balanced 
by the pumping of energy to the stable second over- 
tone that behaves as an energy sink, leading to the full 
limit-cycle stability. 



3 Convective models in 2-D 

3.1 From radiative to convective setups 

After the previous results obtained on purely radia- 
tive models in 1-D, we next address the convection- 
pulsation coupling that is suspected to quench the ra- 
dial oscillations of Cepheids close to the red edge. The 
idea is to slightly modify the unstable radiative setup 
to obtain a convective zone superimposed to the ion- 
isation region in 2-D simulations. The convective in- 
stability obeys Schw arzschild's criterion given by (e.g. 
Chandrasekharl 1 1 96 ll ) 



dTo 



dz 



> 



Ehot > 



gKo{To) 



(5) 



where Fbot = Ko\dTo/dz\ is the imposed bottom ffux. 
For a given gravity g and conductivity hollow Ko{To), 
convection will therefore develop for a large enough bot- 
tom flux satisfying Eq.©. 
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Fig. 5 Snapshot of the vorticity field V xv in i 
with convection and a^- mechanism. 



2-D model 



Figure [5] displays a snapshot of the vorticity field ob- 
tained in such a convective simulation at time t = 2500. 
In order to ensure that the thermal relaxation is well 
achieved, we ran this simulation for a much longer time 
than in the purely-radiative case and the resolution is 
256 X 256 gridpoints. As expected, convection develops 
in the middle of the box where the conductivity hollow 
is locateclll. Moreover, because of the density contrast 
across the convection zone, the vorticity is trapped in 
strong convective plume downdrafts that easily over- 
shoot in the stable zone below due to their low Peclet 
number (^Dintransi 12009. ). 

^An animation of this simulation is provided at 
|http : //www . ast . obs-mip . f r/ users/ tgast ine/ conv_kappa. avi^ 
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3.2 Evolution of acoustic modes with convection 

To address the time evolution of the acoustic modes 
that are present in the DNS with convection, we project 
as before the physical fields onto an appropriate acous- 
tic subspace to get the projection coefficient Qn(^)- 
Here £ denotes the mode degree while n is its order 
(with kx = {27r /Lx)i and is the width of the box). 




(dashed green line). Indeed, after a long transient dur- 
ing which it reaches a non-trifling value (~ 30%), it 
progressively tends to zero at long times while the con- 
vective motions are responsible for the whole kinetic 
energy content. The radial oscillations excited by hz- 
mechanism are thus quenched by convection and this 
situation is relevant to the red edge where the unstable 
acoustic modes are known to be damped by convective 
motions occuring at the star's surface. 

We thus have shown that the Rayleigh number is 
a key parameter for this problem. But further inves- 
tigations of the convection-pulsation coupling are in 
progress, by especially checking time-dependent theo- 
ries of convection for which several free parameters can 
be constrai ned by such kind of direct numeric al simu- 
lations (e.g. Yecko, Kollath, and Buchler 19981 ). 
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Fig. 6 Evolution of the acoustic kinetic energy ratio for 
two simulations with a weak convection (solid blue line) and 
a stronger one (dashed green line). 

Fig. |6] shows the total energy embedded in acoustic 
modes for two simulations with convection that differ 
by their Rayleigh number given by 



Ra = 



ds (V - Vadia)^<^; 
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conv 



(6) 



where Hp = —{dhipQ/dz)~^ is the pressure scale 
height, (iconv the vertical extent of the convection zone, 
X = Kq / pQCp the radiative diffusivity and s the entropy. 
This number allows to quantify the strength of the con- 
vective motions, that is, the critical Rayleigh number 
above which convection d evelops is about 10 ^ for poly- 
tropic stratifications (e.g. ICough et al. 1976 ). 

In the first simulation (Ra = 4 x 10^, solid blue line), 
the Rayleigh number is slightly super-critical and the 
convection intensity is weak. In this case, we see that 
the acoustic energy ratio remains large (i.e. > 80%, 
the remaining 20% being in the convection) and the 
nonlinear saturation is reached in the same way than 
in the previous radiative models. In other words, the 
mode propagation is not affected by convective motions 
and the convection can be seen as frozen-in. 

On the contrary, increasing the Rayleigh number 
leads to a more vigorous convection and the acoustic 
energy ratio is much more affected by convective plumes 
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